/*
 * ScanLargeArrays : Scan is done for each block and the sum of each
 * block is stored in separate array (sumBuffer). SumBuffer is scanned
 * and results are added to every value of next corresponding block to
 * compute the scan of a large array.(not limited to 2*MAX_GROUP_SIZE)
 * Scan uses a balanced tree algorithm. See Belloch, 1990 "Prefix Sums
 * and Their Applications"
 * @param g_out g_out data 
 * @param g_in  g_in data
 * @param block  local memory used in the kernel
 * @param sumBuffer  sum of blocks
 * @param length length of the g_in data
 */

__kernel
void blockAddition(__global int * g_in, __global int * g_out)
{	
	int globalId = get_global_id(0);
	int groupId = get_group_id(0);
	int localId = get_local_id(0);

	__local int value[1];

	/* Only 1 thread of a group will read from global buffer */
	if(localId == 0)
	{
		value[0] = g_in[groupId];
	}
	barrier(CLK_LOCAL_MEM_FENCE);

	g_out[globalId] += value[0];
}


__kernel 
void prefixSum(__global int * g_in, __global int * g_out, __local  int * l_block, const int length)
{
	int thid = get_local_id(0);
	
	int offset = 1;

    /* Cache the computational window in shared memory */
	l_block[2*thid]     = g_in[2*thid];
	l_block[2*thid + 1] = g_in[2*thid + 1];	

    /* build the sum in place up the tree */
	for(int d = length>>1; d > 0; d >>=1)
	{
		barrier(CLK_LOCAL_MEM_FENCE);
		
		if(thid<d)
		{
			int i = offset*(2*thid + 1) - 1;
			int j = offset*(2*thid + 2) - 1;
			
			l_block[j] += l_block[i];
		}
		offset *= 2;
	}

    /* scan back down the tree */

    /* clear the last element */
	if(thid == 0)
	{
		l_block[length - 1] = 0;
	}

    /* traverse down the tree building the scan in the place */
	for(int d = 1; d < length ; d *= 2)
	{
		offset >>=1;
		barrier(CLK_LOCAL_MEM_FENCE);
		
		if(thid < d)
		{
			int i = offset*(2*thid + 1) - 1;
			int j = offset*(2*thid + 2) - 1;
			
			int t = l_block[i];
			l_block[i] = l_block[j];
			l_block[j] += t;
		}
	}
	
	barrier(CLK_LOCAL_MEM_FENCE);

    /*write the results back to global memory */
	g_out[2*thid]     = l_block[2*thid];
	g_out[2*thid + 1] = l_block[2*thid + 1];
}

__kernel 
void ScanLargeArrays(__global int * g_in,
               		 __global int * g_out,
              		 __local  int * l_block,	 	// Size : block_size
					 const int block_size,	 		// size of block				
              		 const int length,	 	 		// elements 
					 __global int * sumBuffer)  	// sum of blocks
			
{
	int thid = get_local_id(0);
	int gid = get_global_id(0);
	int bid = get_group_id(0);
	
	int offset = 1;

    /* Cache the computational window in shared memory */
	l_block[2*thid]     = g_in[2*gid];
	l_block[2*thid + 1] = g_in[2*gid + 1];	

    /* build the sum in place up the tree */
	for(int d = block_size>>1; d > 0; d >>=1)
	{
		barrier(CLK_LOCAL_MEM_FENCE);
		
		if(thid<d)
		{
			int i = offset*(2*thid + 1) - 1;
			int j = offset*(2*thid + 2) - 1;
			
			l_block[j] += l_block[i];
		}
		offset *= 2;
	}

	barrier(CLK_LOCAL_MEM_FENCE);

	int group_id = get_group_id(0);
		
    /* store the value in sum buffer before making it to 0 */ 	
	sumBuffer[bid] = l_block[block_size - 1];

	barrier(CLK_LOCAL_MEM_FENCE | CLK_GLOBAL_MEM_FENCE);
	
    /* scan back down the tree */

    /* clear the last element */
	l_block[block_size - 1] = 0;	

    /* traverse down the tree building the scan in the place */
	for(int d = 1; d < block_size ; d *= 2)
	{
		offset >>=1;
		barrier(CLK_LOCAL_MEM_FENCE);
		
		if(thid < d)
		{
			int i = offset*(2*thid + 1) - 1;
			int j = offset*(2*thid + 2) - 1;
			
			int t = l_block[i];
			l_block[i] = l_block[j];
			l_block[j] += t;
		}
	}
	
	barrier(CLK_LOCAL_MEM_FENCE);	

    /*write the results back to global memory */

	if(group_id == 0)
	{	
		g_out[2*gid]     = l_block[2*thid];
		g_out[2*gid + 1] = l_block[2*thid + 1];
	}
	else
	{
		g_out[2*gid]     = l_block[2*thid];
		g_out[2*gid + 1] = l_block[2*thid + 1];
	}
}

